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This thesis examines the differences of deterministic and stochastic 

LANCHESTER-type combat models. Using an example of square-law attrition, 
solution methods and solutions are described. A new analytic solution 
for equal attrition rate coefficients is given. The numerical comparison 
includes hypotheses about the expected force levels and the variability 
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I. INTRODUCTION 



Combat models are widely used as decision aids in the defense- 
planning process, at least within the NATO alliance. Current operation- 
al combat models are very complex because combat is a very complex 
process. Unfortunately it is difficult (if not impossible) for the 
beginner to understand the modelling approaches, concepts and motivation, 
that may have been used to build such operational models. However, one 
frequently considers a simple model as a paradigm for the development and 
understanding of such complex models. This basic approach will be used 
in this thesis to explore certain issues in the on-going debate about 
the relative merits of stochastic and deterministic combat models. 

A simple model is examined to explore differences between a determin- 
istic and a stochastic approach to a certain type of analytical combat 
model. As already mentioned, combat is a very complex process, but it 
is also a complex random process, which can be supported by many examples 
from military history. Analytical models are abstractions and very often 
simplifications of reality. It seems to be a legitimate question to ask, 
what effects the further abstraction of neglecting the randomness in 
combat may have. At this moment, it should be pointed out that within 
existing operational analytical models, both stochastic and deterministic 
models are used. 

Previous work done by SPRINGALL [9] and CLARK [4] evolved around 
theoretical aspects. Their main concern was to give exact analytical 
solutions and their proofs. CRAIG [5] started to explore the differences 
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between stochastic and deterministic models more from the numerical 
point of view, which will be continued in this thesis. 

In the next chapter, a deterministic and stochastic version of a 
differential combat model will be described. The deterministic versions 
are well-known as LANCHESTER's equations of modern warfare, which were 
developed in 1914. Combat models, which model attrition from enemy 
action through a system of differential equations, are usually referred 
to as LANCHESTER-type models of warfare. 
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II. THE PARADIGM MODELS 



A. THE DETERMINISTIC MODEL 

First, LANCHESTER's equations of modern warfare will be briefly 
reviewed and some simple extensions given. 

In 1914 LANCHESTER [7] hypothesized that under "modern conditions" 
in a combat between two homogeneous forces the firepower of the surviving 
weapons of one side can be concentrated on the surviving targets of the 
other side, so that each side's casualty rate is proportional to the num- 
ber of enemy firers. This can be described by the following equations: 



dx 

dt 


-ay 


(2.1) 


= 

dt 


-bx 


(2.2) 


i initial conditions 




x(o) 


= x 0 


(2.3) 


Y(0) 


= ^0 


(2.4) 



where a is the attrition rate with which the Y-force attrits the 
X-force, similarly for b. Xg and Yg are the initial force levels and 
X(t) and Y(t) are the force levels at time t. The force levels, as a 
function of time t, can be written as 



X(t)=XgCosh(Vab t)-Ya/b ygSinh(Yab’ t) 


(2.5) 


Y(t)=ygCosh(-Vab’ t)-Vb/a XgSinh(^ab t) 


(2.6) 



The state equation relating initial force levels with force levels at 
some time t can be derived by dividing (2.1) by (2.2), separating 
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variables and integrating to yield 

b(x Q 2 -X(t) 2 )=a(y 0 2 -Y(t) 2 ). (2.7) 

This form of the state equation explains why this model is also referred 
to as the "square-law" attrition process. WEISS [11] has given a set of 
assumptions under which LANCHESTER's equations for modern warfare may 
apply: 

A1 . ) Two homogeneous forces are engaged in combat. Every unit 
on a particular side has the same capabilities. The attrition 
rate may be different for the two forces. 

A2.) Each unit on one side is within weapon range of all units 
on the other side. 

A3.) The effects of successive rounds on the target are 
independent. 

A4.) Each unit has perfect knowledge of target locations and 
fires only at live target (one at a time) killing them at a 
constant rate, which does not depend on the number of targets 
alive. 

A5.) Fire is uniformly distributed over surviving targets. 

The above model implies a fight until one force is annihilated. 
Therefore the model will be slightly changed by introducing the concept 
of unit breakpoints, X bp and Y bp , which are simply force levels at which 
the side, who reaches it first "breaks off" the engagement, leaving the 
other side as a winner. Also, to be more precise, it should be noted 
that negative force levels for breakpoints equal zero or force levels 
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less than nonzero breakpoints are impossible. So the deterministic 
LANCHESTER-type combat model with "square-law" attrition takes the 

v x(t) 

y bp <Y(t > < 2 - 8 ) 

otherwise 
*bp <x(t) 

y bp <v(t) (2.9) 

otherwise 
with initial conditions 

X(0) = x Q (2.10) 

V(0) = y Q (2.11) 

The model in this form, equations (2.8) through (2.11), was used for 
comparisons throughout the thesis. 

B. THE STOCHASTIC MODEL 

There are several ways to include random variations in LANCHESTER- 
type models. These are: 

* The attrition rate coefficients may be random variables. 

* The enemy's initial force level may be a random variable, 
weakening the assumption of perfect knowledge. 

* The breakpoints may be random variables. 

* The casualty rate is fixed, but the occurrence of casualties 
over time may be random. 

The only random variation considered here will be the random occur- 
rence of casualties over time. Another specification was to choose a 



following form: 

-ay 



dx 

at 



= < 



p. = 

dt 



0 

-bx 

0 
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model similar to the "square-1 aw"attrition in order to allow comparisons. 

In other words, the question to be asked is "how do random fluctuations 
in the occurrence of casualties modify the deterministic results of the 
square-law attrition process?" 

The approach used here was a continuous parameter MARKOV chain model, 
where the time t varies continuously and the number of combatants on 
each side is a non-negative integer. Let M(t) be the size of the X-force 
at time t with a particular state value m. Let N(t) be the size of the 
Y-force at time t with a particular state value n. Let iHq and ng be the 
initial force levels and m^, n bp be the breakpoint force levels of X 
and Y respectively. Fig. 1 shows the state space of this MARKOV chain 
model. Note that at a given time t, any state is described by the two 
force levels of the X and Y force. As each side loses units due to 
attrition and no replacements are allowed, it is easy to understand why 
BILLARD [1] referred to this type of process as a "bivariate death process. 
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Figure 1 - STATE SPACE OF THE MARKOV CHAIN MODEL 
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For the description of the stochastic square-law attrition process 
corresponding to the two deterministic differential equations (2.8) and 
(2.9), a system of many differential equations, depending on the battle 
termination model, is required. This system will be given for a fixed- 
force-level-breakpoint battle with square-law attrition. 

The following assumptions yield the stochastic square-law attrition 
process. 

Al.) The attrition process depends only on the current system 

state and time and not on the past history (this assumption is 

usually referred to as markovian property). 

a-? ^ pf one x casualty during the \ _ - 

A2.) The probability P^ tjnle 1nterval , t t0 t ? h ) - ah 

A3 1 The Drobabilitv P ( one y casualty during the ) m .. 

m.) me proDaoi nty i nter val t to t+h / bh 

A4.) The probability of more than one casualty occurring in 
the time interval t to t+h is of the order of magnitude o(h), 
where lim^Qo(h)/h = 0. 



A5.) No more casualties can occur once m = m^ or n = n^. 

Making the time interval h infinitesimally small, the following 
set of forward CHAPMAN- KOLMOGOROV equations can be developed. Let 
P(t,m,n) be the probability that the system is in state (m,n) at a 
time t. For convenience each equation is related to a region in the 
state space shown in Fig. 2. 
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( 2 . 12 ) 



For m=mg and n=n Q , Region I 
d P 

j^(t,m 0J n 0 ) = -(an 0 +bm 0 ) P(t,m 0 ,n 0 ) 
for m b p<m<mQ and n = ng , Region II 



d P 

^■(t,m,n 0 ) = an Q P ( t ,m+l ,n Q ) - ( an Q +bm) 
for n bp <n<n 0 and m=m 0 , Region III 


P(t,m,n Q ) (2.13) 


d P 

j-£-( t ,m Q ,n ) = bm Q P ( t ,m Q ,n + l ) - ( an+bm Q ) 
for m t)p <m<ni 0 anc * n bp <n<n 0’ Re 9 lon IV 


P(t,m 0 n) (2.14) 


d P 

j-jr(t,m,n) = anP ( t ,m+l ,n )+bmP ( t ,m ,n + l ) ■ 
for and n bp <n<n 0’ Pe 9 10n V 


• ( an+bm ) P ( t ,m ,n ) 

(2.15) 


d P 

Jjr(t, m b p>n) = anP(t,m bp + l ,n) 
for n=n b p and m 5 p <m<m Q> Region VI 


(2.16) 


d P 

dt( t,m,n bp) = bmP (t,m,n bp +l ) 
for m=m b p and n=n bp , Region VII 


(2.17) 


P(t,m bp , n bp) = 0 for all t 


(2.18) 


because of the definition of a breakpoint 
initial condition is 


force level. The 


P ( 0 ,mQ ,ng ) - 1.0 


(2.19) 
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Figure 2 



REGIONS IN THE STATE SPACE 
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As P(t,m,n) is a joint probability distribution, the follow- 
ing must also be true 



0- P ( t ,m ,n ) - 1 . 0 



and 



m. 



2 H 

m=m bp n=n b P 



P(t,m,n) = 1 .0 



( 2 . 20 ) 



( 2 . 21 ) 
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III. SOLVING THE DIFFERENTIAL EQUATIONS 



A. SOLUTIONS FOR THE DETERMINISTIC MODEL 

Force levels as a function of time were already given in equations 
(2.5) and (2.6). It is relatively easy to obtain analytical solutions 
for such simple deterministic models as described before. On the other 
hand, it is necessary to point out that for models with any degree of 
operational realism, analytical methods for solving the differential 
equations are usually not available. Therefore, some numerical method 
with a digital computer is usually used. TAYLOR [10] has summarized in 
his Appendix C the most widely used numerical methods, a discussion of 
which seems unnecessary at this point. 

B. SOLUTIONS TO THE STOCHASTIC MODEL 

Even for this relatively simple stochastic model with fixed force 
level breakpoints, which are usually nonzero, a complete set of general 
analytical solutions for the differential equations (2.12) through (2.19) 
has not been found. On the other side, given some minor restrictions 
like a fight to the finish or equal attrition rate coefficients, solu- 
tions, or at least solution methods have been proposed which will be 
briefly discussed in the next section. 

First, the method for getting the state probabilities used here will 
be described. Numerical solutions were obtained using the fourth-order 
RUNGE-KUTTA method, which is probably one of the best known finite dif- 
ference approximations to ordinary differential equations (next to the 
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EULER-CAUCHY-method) . To increase the accuracy of the overall solutions 
analytical results for certain regions of the state space were substituted. 
These analytical solutions will be stated now. For region I, i.e. no 
casualties on either side, the solution to (2.12) is 

P(t,m 0 ,n 0 ) = exp -(an 0 +bm 0 )t (3.1) 



which can be derived by the standard method used for this kind of dif- 
ferential equation. For the boundary cases, region II and III, i.e. 
one of the two sides has not yet had a casualty, TAYLOR [10] has given 
the analytical expression as 



for m 5p <m<rn Q and n = n n , Region II 



P(t ,m. 



n 0 } J! | ; 



'0 



an 0 /b(e bt -l) 



, exp(-(bm 0 +an 0 ) t) , (3.2) 



where J = trig-m 



for n k <n < n g and m=mg, Region III 



P(t 



» m 0 » n ) " K! j 



F, {bm Q / a ( e at - 1 ) 



, exp(-(bm 0 +an 0 ) t) , (3.3) 



where K = ng-n. 

These two equations (3.2) and (3.3) were obtained by 
recursively solving equations (2.12), (2.13) and (2.14) "from 
the top down." 



C. OTHER ANALYTICAL SOLUTIONS 

The solutions or solution methods for getting the state 
probabilities will only be stated for the square-law attri- 
tion process. Only two were used for the numerical work for 
this thesis. 
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Apparently, one of the "oldest" analytic solutions was 



given by BROWN [3] in 1955 for the general stochastic 
LANCHESTER- ty pe combat model with time independent attrition 
rates. His approach and solution will be briefly outlined 
for square-law attrition. Consider a path from state 
( m Q , n Q ) to some state (m,n). This path can be described 
as a sequence of J = ( nig -m) zeros and K = ( ng -n) ones, 
where a zero corresponds to an X casualty and a one to a Y 
casualty. Using the binary representation of a positive 
integer, one can define to each realization of a battle path 
an integer k given by 



corresponding to k is a Y casualty and d. = 0 otherwise. 

K j i 

Also let Ij K be the set of all positive integers whose 



binary representation contains exactly K ones and J zeros. 
Then 



k d k , 1 d k ,2 ' * ' d k , J + K 

L. 

where d. = 1 if the r casualty along a battle path 




r 



k , r ‘ n 0 




Then BROWN [3] has shown that 
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du , (3.4) 



1 «? J + K- 1 

P(t,m,n) ”27 / PC K k ,r 

keI J.K •" r=0 
exp(-iut)-exp(-tl(m kjr ,n kjr )) 

1 ^ m k,r’ n k,r^“* u 

where i = V-T , 

1 (m ,n ) = an+bm , 

K _ g k ,r+1 

k,r 1 - i u / 1 ( m k>r > n k >r ) 

and 

9 k , r+1 = d k,r+l an k,r +(1 “ d k,r+l )bm k,r * 

There was no indication that nonzero breakpoints were 
excluded. A discussion of this solution follows in the next 
section. 

About 14 years later, in 1969, CLARK [4] proposed another 
approach which TAYLOR [10] called a "hybrid analyti cal -numer- 
ical" method. The restriction is that the breakpoints have 
to be zero, i.e. it is a fight to the finish. Although 
proposed for a general time independent attrition function, 
this approach will be outlined for "square-law" attrition. 
Then according to CLARK [4] the state probabilities are given 
by 
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for 0<m<m 0 and 0<nln Q 



m ° ” Q 

P(t,m,n) = 2 _ c j’k • exp ^-(ak + bj ) tj , (3.5) 



j=m k=n 



for 0<m<m 0 and 



n = 0 



u u_ 

P( t ,m ,0) = cj!’° + H ZL C*’° • exp (- ( ak + bj ) t) , (3.6) 
u,u j-m k=l J ’ K v ' 



for m=0 and 0<n<n 



0 



P ( t , 0 , n ) = + 



m 0 


"0 


E 


Z 


j=m 


II 


m 0 


n 0 


II 


j=m 


II 

-V 



C j]k * ex P (- ( a k+b j ) t) ,(3.7) 



and at last thereis to remember that P(t,0,0) = 0 for all 
times t. The constants C^*[J are determined by a system of 

J i K 

partial difference equations. 



For 0<m<j<mQ and 0<n<k<n 



0 



r m , n 
j >k 



r m+l ,n , r m ,n + l 

an L j ,k bm ,k 

a(n-k)+b(m-j) 



, (3.8) 
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for 0<m<j<niQ and 0<n=k<n Q 




J .n 



an 



m+1 ,n 
C . 

J 



b(m-j ) 



for 0<m=j<niQ and 0<n<k<ng 



m,n m,n+l 

C . = bm C 

m j k m j k 

a(n-k) 



for 0<m=j<mg and 0<n=k<ng 
but (m,n) f ( m 0 , n 0 ) 

m,n m 0 m 0 m,n ^0 - m,n 

C = - yz y_ c. . “ / C . . 

m ’ n j=m k = n+l J,k j =m+l J ’ k 

with c!°’"° = 1.0. 

m o ,n 0 



, ( 3 . 9 ) 



, ( 3 . 10 ) 



, ( 3 . 11 ) 
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Also for 0<m<j<mg and 0=n<k£ng 



r m,0 

L j,k 



bm 



, 1 
j,k 



ak+b j 



similarly for 0=m<j<m 0 and 0<n<k<n Q 



(3.12) 



r 0 ,n 
L j,k 



p 1 > n 
anC .i,lc 

ak+b j 



(3.13) 



Then for l<mlmg 



».0 >- Jill. m>0 

*t>.0 ' -2_ /— C j.k 

j=m k=l 



(3.14) 



and finally for lining 




j=l k=n 



Though having the publishing date of 1979, the next 
approach was published in June 1980 by BILLARD [1]. She 
considered the LANCHESTER- type square-law attrition combat 
model as a pure death-process and applied SEVERO’s [8] 
recursive theorem for solving differential equations. As 
before, only a fight to the finish has been considered. 
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The first step is to identify each point (m,n) in the 
state space by a counting coordinate k, where 
k = ( m 0 +l)( n Q +l) 

-m( n Q +l ) -n . (3.16) 

Then 

P [t ,m ,n] = P[t,k] (3.17) 

and the differential equations (2.12) through (2.18) take on 
a slightly different form. As an example, (2.15) will be 
given by 

H P 

j£-(t,k) = anP [t ,k-n Q -l ] +bmP [t ,k-1 ] 

- ( an+bm) P [t , k] . (3.18) 

The whole set of differential equations was then 
expressed in matrix terms as 



^•P(t) = BP(t) 


(3.19) 


with a solution given as 




P(t) = Ce (t) 


(3.20) 



where ejt) is the (mg+1 ) (hq+ 1 )xl -vector with elements 

j. L. 

exp(bj c t) with b^ being the k tn diagonal element of the 
matrix B^. The matrix B_ can be partitioned into submatrices, 
whose m-coordinate is common, due to the ordering defined by 
the counting coordinate k (Equation 3.16). 

Then 

B = (b uv ) , u ,v = 1,2.. . m Q +l 
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where the submatrices b uv have the elements 
-uv = ( b uv ^ p ’ q ^’ p,q = 1 »2...n 0 + l. 

So the matrix has the elements 
b uu (P>P) = -a(n Q -p+l )-b(m Q -u + l ) 

for u = 1,2. ..mg and p = l,2...n 0 

b uu (p,p-l) = b(m 0 -u+l) for u = l,2...m 0 +l 

p = 2,3... n n +l 

and 0 

b u> u-l( p ’ p ) = a(n Q -p+l) for u = 2,3...m Q +l 

P = 1,2... n Q +l . 

All other elements are zero. 

Thus, the matrix B_ has at most 3 nonzero entries per row 
or column. The matrix C_ can be partitioned in the same way. 
Then using SEVERO's [8] theorem and the special form of the 
matrix B^, only a part of the £-matrix needs to be determined. 
This part will be omitted here, but the final result will 
be given by 

k 

P ( t , k ) = y ' c ( k , j ) . exp ( b j t ) (3.21) 

j = l 

x. L. 

where c(k,j) is the (k,j) L element of the solution matrix £. 

The previous two approaches have required that the force 
level breakpoints be zero. Now, a result will be given whose 
restriction is that the attrition rate coefficients be equal, 
but nonzero breakpoints are allowed. For further reference 
it will be called the Equal -Attri ti on-Ra te-Coef f i ci ent- 
Solution (EARCS). 
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Let a = b = f. 



For m bp <ro<m 0 and n bp <n<n 0 



Pfr m nl = c(m.n) 

’ ’ (m 0 + n o-tn- n )! 



(i. e - ft ) ra o +n o- m - n 



exp(-f (m+n)t) 



(3.22) 



for m = n bp and n bp <n<n Q 



P(t.m bp .n) = fn 

J ! 




l-exp(-ft(m b +l+n+k)) 



k = 0 



bP 

f ( m bp + l +n + l< ) 



')■ 



(2.23) 



where J = m o +n o" m bp" 1 _n 
for n = n bp and m bp <mFm 0 



P(t,m,n bp ) = fm c ^ m>n bp +1 l 

K ! 



(-D 



j=0 



*(;) 



1 -exp (-ft(m+l+n bp +j )) 



f (m+1 +n bp +j ) 






(3.24) 



where K = m Q +n 0" m " n bp” 1 • 

The coefficients C(m,n) satisfy for m b p<m<mQ and n bp <n<n Q 
the following partial difference equation 

C(m,n) = nC(m+l,n) + mC(m,n+l) (3.25) 
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with the boundary conditions 

m^-m 

C(m,n 0 ) = (n Q ) for m bp <m<m 0 



and 



n 0“ n 

C(m Q ,n) = (m Q ) for n bp <n<n Q 

This result has been developed using a method verbally pro- 
posed by TAYLOR. The method will now be outlined. ISBELL 
and MARLOW [6] described a stochastic LANCHESTER-type attri- 
tion process with a different attrition function. Instead of 
the attrition of one force being proportional to the number of 
of enemies of the opposing force as in the square-law attri- 
tion (e.g. for the attrition of the M-force let 
A(m ,n ) = an 

be the attrition rate and 
B(m,n) = bm 

the attrition rate for the N-force with square-law attrition), 
their attrition rates looked like 
A(m,n) = an+cm 



and 



B(m,n) = bm+dn 
with the restriction that 
a+c = b+d. 

But with c=d=0 and a=b=f we are back to square-law attrition. 
This leads to equations (3.22) and (3.25). Equations (3.23) 
and (3.24) were derived in the following way (e.g. for 3.23). 
Solving equation (3.22) for m = m bp +l yields 
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(3.26) 




0 "0 "'bp 



. (l_e _ft ) n, 0 + n 0" m bp‘ 1 ‘ n 



substituting for the second factor its BINOMIAL expansion 



with J = m o + n Q -rn 5 p“ 1 an d multiplying through by the third 
factor. Then using equation (2.16), the differential equation 



into the extended form, it can now easily be integrated to 
yield equation (3.23). 

D. DISCUSSION 

In the discussion of the analytical solutions outlined in 
the last section, there is one important point. BROWN [3] 
himself points out that unless m is close to mQ and n is 
close to ng his result (equation (3.4)) is of "little practi- 
cal interest." Most of the analytical sol uti ons , especi al 1 y for 
more general LANCHESTER- type models, have little more than 
"symbolic" character. BROWN's solution is a good example 
of that. 

In comparing the solutions given by BILLARD [2] and 
CLARK [4], this author has the feeling that both solutions are 
equivalent and only the representation is different. This 




k=0 



e- ftk (l) J ' k 



for m = m^p and n bp <n - n g , and substituting equation (3.26) 
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intuitive guess needs verification. It may be concidence 
that CLARK [4] and SEVERO [8] published their work in the 
same year. 

The last presented solution (equation (3.22) through 
(3.24)) seems to be relatively handy for use on a computer. 

It has a big advantage over numerical solution methods other 
than its accuracy, because it is an exact result. Like 
CLARK's method, the coefficients have to be calculated only 
once for a given set of input data. Then, to get the state 
probabilities for a certain point in time you have to make 
only one set of calculations, as opposed to the numerical 
methods where one has to go from time t = 0 to time t = t 
in small time steps and then have only an approximate result. 

E. IMPORTANCE OF THE STATE PROBABILITIES 

The state probabilities as a function of time are the 
key to calculating several quantities of interest. These 
are expected force levels as a function of time, variances 
and standard deviations in the force levels and also the 
probability of winning. These quantities are necessary to 
legitimately compare the stochastic with the deterministic 
resul ts . 

To get at least a feeling of how the state probabilities 
evolve over time, the joint probability distribution will be 
presented in a 3-D-picture. It is indeed surprising that 
more use has not been made of computer graphics to 



29 



investigate the dynamics of a stochastic LANCHESTER- type 
combat model. Table 1 gives the data used for the next five 
figures. These figures may be thought of as "snapshots" of 
the joint probability for the survivors in this battle taken 
at a sequence of increasing times. 
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TABLE 1 



Data for the Numerical Example 1 



Force Levels 


m o 


= 40 


m. n = 0 
bp 






n 0 


= 40 


Hur. = 0 
bp 




Attrition Rates 


a = 


0.008 


M casualties 
minute and N 


per 
f i rer 




b = 


0.004 


N casualties 
minute and M 


per 
f i rer 



At Times t-j = 0.025 tf 

t 2 = 0.250 t f 
t^ = 0.500 tf 
t^ = 0.750 tf 
t 5 = 1 .000 t f 



where tf = 155.81 minutes is the time a 
deterministic battle with the same force 
levels and unit breakpoints ends, i.e. 



x(t f ) = 0.0 and y(tf) = 28.28 
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'k P (t,m,n) 




Figure 3 - PLOT OF JOINT PROBABILITIES P(t,«,n; 
with data according to Table 1 
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H p(t,m,n) 




figure 4 - PLOT OF JOINT PROBABILITIES P(t,«,n / 
with data according to Table 1 
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1 > P(t ,m,n) 




Figur* 5 - -PLOT OF JOINT PROBABILITIES P(t,»,n; 
with data according to Tablw 1 
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ii P(t,m,n) 



Force Level, n 




Force Level, n 



figure 6 - PLOT OF JOINT PROBABILITIES P(t,«,n) 
with data according to Table 1 
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, i P(t,m,n) 




rigur* 7 - PLOT OF JOINT PROBABILITIES P(t,«,n) 
vith data according to Tabl# 1 
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At the time t = 0 all probability is located at 
(itiq.iIq) in the state space (Region I). As time passes, the 
probability mass is distributed over more states, with the 
mode moving away from the diagonal towards the winning side. 
All points in the state space with breakpoints, i.e. (m b p,n) 
and (m,n b p) for all m and n, are absorbing states, probability 
mass is absorbed in that states. The sum of probability mass 
in Region V (see Fig. 2) represents the probability, that the 
N-force wins at that given time, in Region VI that the M-force 
wins. 

The next five figures show a similar sequence of plots for 
the joint probability P(t,m,n) for the force levels M(t) and 
N(t). The data is explained in Table 2. Note the small dif- 
ferences because of the nonzero force level breakpoints. 
Probability mass having reached the breakpoint "piles" up 
there. The state space is reduced by the fixed force level 
breakpoints, but the probability distribution evolves basical- 
ly in the same qualitative manner as in the previous example. 



37 



TABLE 2 



Data for the Numerical Example 2 



Force Levels 



m 0 = 40 m bp = 8 

n 0 * 40 "bp * 8 



Attrition Rates a = 0.008 M casualties 

minute and N 

b = 0.004 N casualties 
minute and M 



At Times t-j = 0.025 t^ 

t 2 = 0.250 t f 
t^ = 0.500 t^ 
t 4 = 0.750 t f 
t 5 = 1 .000 t f 



where t^ = 120.68 minutes is the time a 
deterministic battle with the same force 
levels and unit breakpoints ends, i.e. 
x(t^) = 8.0 and y(t^) = 28.83 



per 
f i rer 

per 
f i rer 
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^ P(t,m,n) 




Figura 9 - PLOT OF JOINT PROBABILITIES P(t,a,n) 
with data according to xabl# 2 
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J> P(t,m,n) 




Force Level, n 



Force Level, m 



figure 1u - PLOT OF JOINT PROBABILITIES P(t,a,n; 
with data according to iable 2 
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i k P{t ,m,n) 



Force Level , m 




Force Level , n 



Figure 11 - PLOT OF JOINT PROBABILITIES P(t,«,n; 
with data according to Table 2 
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P(t,a,n) 




jrigura 12 - PLOT OF JOINT PROBABILITIES P(t,«,a) 
with data according to Tabla 2 
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IV. EVOLUTION OF THE FORCE LEVELS 



A. THE DETERMINISTIC MODEL 

The evolution of force levels as a function of time, X(t) and Y(t), 
has already been stated in equations (2.5) and (2.6). The next two 
figures show the force levels with different breakpoints. It is easy to 
realize that introducing a nonzero breakpoint does not change the under- 
lying function. In other words, using X(t) as an example, for both 
figures the same curve was used but at the point where the X-force reaches 
its breakpoint, the curve is "cut." From that point in time, there are 
no more changes in the force levels. So introducing a nonzero break- 
point only shifts the discontinuity (marked by DX1 in Fig. 13) up along 
the curve to DX2. 

The probability for one side to win is either one or zero, because 
it is a deterministic model. To easily determine which side is going to 
win, a victory prediction condition can be obtained by solving each force 
level equation (2.5) and (2.6) for the time to reach its breakpoint tx bp 
by substituting X(t)=x bp and ty bp by substituting Y(t)=y bp . Then X will 
win if ty bp < tx bp , which leads to the prediction condition. X will win 
a fixed force breakpoint battle if and only if 

. 

X ° y0> bUo-Xbp) • < 4 -’> 

This shows that given the initial data one can predict the outcome 
of the battle in terms of force levels and time until the battle 
finishes . 
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Figure 1 3 




- FORCE LEVELS OVER TIME 

Deteraini6tic Model with Breakpoint x ^p = 0 
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X, 



Figure 14 




- f'ORCE LEVELS OVER TIME 



Deterministic Model with Breakpoint 



*b P > 0 
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B. THE STOCHASTIC MODEL 

Every possible state (m,n) in the state space has associated with it 
a certain probability between zero and one which is a function of time. 

In order to gain more insight into the stochastic process not a single 
realization of a battle has to be considered, but an average battle. 
Therefore the expected value of the force levels (i.e. averages) as a 
function of time and the variances in the force levels were investigated. 

The straight forward way to compute the expected force levels involves 
the knowledge of the probability distribution P(t,m,n). Then 

m 0 

E(M(t) )=y~ m 

m bp n bp 



1 



P(t,m,n) 



(4.2) 



n 0 m 0 

E(N(t) )= P(t,m,n) (4.3) 

n bp m bp 

are the expected values of the force levels as a function of time. 

There are some other ways to calculate the expected force levels, 
one of which will be stated here. Recall the "hybrid-analytical- 
numerical" method to get the state probabilities (equations 3.5 through 
3.15). CLARK [4] has also shown that the i c moment of, for example, 
the M-force level may be computed as 

E(M i (t) ) = d! 1 ! +ZI TL d! 1 ^ exp(-(ak+bj)t ) (4.4) 

u,u j=l k=l J,K 



with for 1-j-mQ and Uk^ 
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(4.5) 



)(!) = 
j ,k 



m 



m=l 



n=0 



r m 9 n 
j ,k 







r 0,n 

L 0,0 



(4.6) 



This again emphasizes the strong point of CLARK's solution. For a 
given set of battle parameters the coefficients C^’[) have to be computed 
only once. Then with this information and relatively small computa- 
tional effort not only the state probabilities but also the first and 
second moment of the force levels can be computed. This determines the 
variance in the force levels at the same time, e.g. 

Var(M(t) ) = E(M 2 (t) ) - E(M(t)) • E(M(t) ) . (4.7) 

On the other side, the weak point is that CLARK (4] considered only 
breakpoints equal zero. 

Many authors have discussed one side's probability of winning or 
probability of winning conditioned on a certain number of survivors, 
which always eliminated the parameter time by integrating from time t=0 
to infinity. This may be legitimate to answer absolute (meaning time 
independent) questions about who will win, but for direct comparisons 
with the deterministic model, this author has the feeling that the best 
question to ask regarding a winner is: 

What is the probability of one side winning given the stochastic battle 
lasted as long as the deterministic one? 

The calculation of these probabilities gives another interesting 
probability, because given the time t = t^. (final time of the determin- 
istic battle) 
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P(battle has not yet finished) = 



1 -P (M wins | t=t^)-P(N wins|t=t^) 



(4.8) 



There has also been work done regarding the distribution of the 
time to finish a battle. But this is beyond the scope of this thesis 
(SPRINGALL [9]). 

C. DIFFERENCES IN THE FORCE LEVELS 

The deterministic model, especially the force level equations (2.5) 

and (2.6), describe a process with continuous state parameters where, in 

reality, the possible states are integers. Quoting from LANCHESTER [7]: 

Since the forces actually consist of a finite number of finite units 
(instead of an infinite number of infinitesimal units) the end of the 
curve must show discontinuity, and break off abruptly when the last 
man is reached; the law based on averages evidently does not hold 
rigidly when the numbers become small. 

LANCHESTER suggested that his differential equations may be good 
approximations only as long as the force sizes are large. He also stated 
that the equations are based on averages, implying an underlying stochas- 
tic process. 

This shows that there must be a difference in the force levels which 
should become significant when the number of combatants is small. This 
difference was called bias by CLARK [4] and TAYLOR [10]. It can be shown 
that 



= -aE(N(t))+aB n (t) 



(4.9) 



and 



dE _ (N( _ t . )j _ s _ bE (M( t ))+ttB m (t) 



(4.10) 



where 
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(4.11) 



m o "o 

B n (t) ’ "bp ZZ P(t ’ ra ' n bp )+ ^Z nP(t ’ m bp" n) 



V' 



V 1 



and 



■0 m o 

B m (t) = >2Z p(t ’ m bp* n »XZ mP(t,m,n b ) . (4.12) 



ii, +1 
bp 



m. +1 

bp 



The bias terms B m (t) and B n (t) can be interpreted as the expected 
values of M(t) or N(t) conditioned on the fact that the battle has 
already ended at time t, for example, 

B m (t)=E [M(t))|N(t)=n bp or M(t)=m bp j (4.13) 

In other words, equation (4.9) says the expected casualty rate of the 
M- force is proportional to the expected number of survivors of the 
N-force given neither of the two forces has reached its breakpoint. 

Define the bias of the X-force as Ax(t) = E(M(t))-X(t) and the bias 
of the Y-force as Ay(t) = E(N(t) ) - Y ( t) . Then using equations (2.8) and 
(2.9) together with (4.9) and (4.10) it follows that 

ax = -aAy+aB n (t) (4.14) 

Ay = -bAx+bB m (t) (4.15) 

with the initial conditions x(0)=) and Y(0)=). This has the solution 

t 

Ax(t)=Vab J {B n (s) •/a/bcosh( Vab* (t-s) 

-B m (s)sinh( Vab (t-s))}ds, (4.16) 
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and 



L _____ 

Ay(t)="/ab J (B m (s) Vb/acosh(-\/ab(t-s) 

-B (s)sinh( Vab*(t-s))}ds 



(4.17) 



Since for a fixed, nonegative argument z, the cosh(z) is always 
greater than the sinh(z), it is easy to visualize that in most of the 
cases both biases are positive, meaning the expected force levels of the 
stochastic model are higher than the deterministic force levels. This 
has been shown by CLARK [4] and CRAIG [5] and confirmed by this author. 
In the rest of the cases the winner's bias is negative or close to zero 
and the loser's bias is positive. Two examples are given in Table 3. 
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TABLE 3 



Examples for cases where one bias is positive and the other 
bias is negative. 

A. Y-force wins in a fight to the finish 



X 

o 

II 

<< 

o 

II 

-P» 

o 


a = 0.08 


x bp = 0 
*bp = 0 


b = 0.04 


t f = 15.581 


Ax ( t^) = 3.22 




Ay(t f ) = -0.23 



B. Y-force wins in a fight with equal initial force levels, 
but different breakpoints. 



X 

o 

II 

<< 

o 

II 

CJ1 


a = 0.08 


x bp - 12 
^bp - 6 


b = 0.08 


t f = 2.85 


ax ( t.p) = 0.15 




Ay(t f ) = -o.io 
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An interesting point has to be mentioned regarding case B of 
Table 3. Here the expected force level of the winner is 
smaller than the expected force level of the loser at the 
time a deterministic battle ends. 

To get a better feeling for the differences in the force 
levels, the next two figures show as an example a large spec- 
trum of force level behavior. The data and notation is 
described in Table 4. There are four battles with four dif- 
ferent breakpoints drawn as they evolve over time from zero 
to the time an equivalent deterministic battle ends. 



i 
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TABLE 4 



o 

II 

O 

X 


x b P (i> 


= 40-0. 2i 


o 

II 

o 

>- 


Y bp<’) 


= 40-0. 2i 


a =0.09 


for i = 


1,2, 3, 4 


b =0.07 






Mi = E M(t) 
Ni = E N(t) 


for i = 


1 ,2,3,4 



Where the index i corresponds to the battle with 

J. L_ 

the i c breakpoint force level: 

X-force is always loser 
Y-force is always winner 
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Figure 1i> - DETERMINISTIC AND EXPECTED FORCE LEVELS FOR 
THE X - FORCE 
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Figure 16 - DETERMINISTIC AND EXPECTED FORCE LEVELS FOR 
THE Y - fORCE 
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CRAIG [5] has formulated hypotheses concerning the biases in the average 
force levels based on his work. His hypotheses were partially confirmed, 
but in some cases they have to be modified. Therefore, another similar 
set of hypotheses will be given and supported by Fig. 15 and Fig. 16, as 
well as some of the later figures. 

H 1) Given fixed initial force levels and attrition rate coefficients, 
as the breakpoint force levels increase, the numerical bias for the 
loser decreases. The biases for the winner do not show this monotone 
trend except for the case of symmetric parity. 

H 2) Everything else constant, the bias of the loser increases with 
increasing initial force levels; this is also true in the symmetric 
parity case for both forces. 

H 3) Given the initial force level ratio is close to one at the time 
corresponding to the end of the deterministic battle, the bias of the 
loser is always larger than the bias of the winner. 

H 4) At the time corresponding to the end of a deterministic battle, 
the biases become larger as the forces come closer to parity. 

H5) The biases at times corresponding to less than one half the 
duration of the deterministic battle are negligible. 

The case of symmetric parity, i.e. equal initial force levels, 
breakpoints and equal attrition rate coefficients seems to be kind of a 
"limiting" case. For example, at parity the biases of both forces behave 
in the same manner and are equal. In Fig. 17 the biases at the time a 
deterministic battle ends as a function of the initial force levels and 
as a function of the breakpoints are presented. It is also another 
verification for the hypotheses H 1 and H 2. 
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Bias 




Flgur# 17 - BIASES IN SYMMETRIC PARITY 
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Figur* 18 - BIASES WITH DIFFERENT INITIAL FORCE LEVELS 
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5.0 
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4.0 4 



3,0 



2.0 - 



l.o - 




x o=yo c 40 
X 8P ® ^BP 

Case I a?b = O.os 
Case Eq^O.oS 
b*0.O7 

CaseHq s o.o8 
b= o.o^ 



0.z 0.4 0.b 0.8 



~ x »7, 



Xo 



Fi«ur« 19 - BIASES WITH DIFFERENT ATTRITION RATE 
COEFFICIENTS 
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Considering the changes in magnitude of the biases for battles like 
in Fig. 15, 16 and 17, CRAIG [5] came to the conclusion that when the 
forces are closer to parity, the biases at the deterministic battle's 
end increase. Several battles were fought starting with symmetric parity 
and then varying the force levels and the attrition rate coefficients. 
Sample results are shown in Fig. 18 and Fig. 19. In Fig. 18, the initial 
force levels were changed giving the Y-force ten and twenty percent 
higher initial force levels. The biases at the end of the deterministic 
battle are plotted as a function of the breakpoint force level ratio 

r = V x o- 

Fig. 18 supports CRAIG 1 s hypothesis, as does Fig. 19. Here not the 
initial force levels but the attrition rate coefficients were changed in 
order to deviate from symmetric parity. The last way to deviate from 
symmetric parity is a case where the deterministic model gives equal 
answers for different battles. The data and the results are shown in 
Table 5. 
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TABLE 5 



Equal initial force level battle with non equal 
X Q = 15 Y q = 15 a = b = f = 



x bp 


Y. 

bp 


AX ( tf 


12 


12 


1 .08 


12 


9 


0.71 


12 


6 


0.15 


12 


3 


0.68 


9 


9 


1.62 


9 


6 


1.14 


9 


3 


1 .06 


6 


6 


2.15 


6 


3 


1 .67 


3 


3 


2.84 



breakpoints 

0.08 

Ay ( t f ) 



1 .08 
0.63 
- 0.10 
0.58 
1 .62 
0.93 
0.75 
2.15 
1 .22 
2.84 
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The X-force level at time t f is the same as the Y-force level. Thus, 
the expected force level for the winner (i.e. the Y-force) is smaller than 
the expected force level of the loser (the X-force), because Ax(t f ) is 
always larger than Ay(t^). This is easy to see when one remembers the 
way the expected force level is computed (equation (4.3)). Since the 
Y-force has a lower breakpoint, there are states (m,n) possible where 
n bp <n<m bp* Apparently these states have a nonzero probability associated 
with them, which decreases the expected force level below the expected 
value for the X-force. This might be a starting point for further 
studies. 

D. VARIABILITY IN THE FORCE LEVELS 

Naturally in the deterministic case there does not exist any variabil- 
ity in the force levels. On the other hand, for the stochastic model, the 
variance in the force levels as a function of time is a measure of the 
dispersion of the number of survivors about their mean value. 

CLARK [4] has hypothesized two different types of behavior for the 
variance in the force levels, which are shown in Fig. 20 for the data 
presented in Table 2. The first type of behavior is that of the variance 
for the N-force, i.e. the variance increases monotonely as a function of 
time and is asymptotic to a limiting value. It was found that this type 
of behavior occurs when the side is going to win and for the case of 
symmetric parity. The second type of behavior shown is the variance of 
the M-force, Var(M(t)), as a function of time. This increases to a 
maximum value then decreases asymptotically to a limiting value. This 
type of behavior is associated with the loser of the battle. 
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The variance in the force levels is a function of the initial and 



breakpoint force levels, the battle time and the attrition rate 
coefficients. Unfortunately one does not know what this dependence is. 
Based on many numerical results a set of hypotheses will be stated and 
the next figures will support them. 

H 1) Given fixed initial force levels and attrition rate coefficients, 
as the breakpoint force levels increase the variance of the force levels 
decrease. 

H 2) Everything else constant the variance of the loser's force level 
increases with increasing initial force levels. This is also true for 
both variances in the case of symmetric parity. 

H 3) Given the initial force ratio is close to one at the time 
corresponding to the end of the deterministic battle the variance 
of the loser's force level is always smaller than the winner's variance. 

H 4) At the time corresponding to the end of the deterministic battle 
the variance in the loser's force level increases as the forces come 
closer to parity. This trend is not true for the variance of the win- 
ner, except for the case of symmetric parity. 

It was the intention of this author that the set of cases to illus- 
trate the hypotheses are the same as in the illustrations (Fig. 17, 18 
and 19) of the hypotheses about biases. So, Fig. 21 shows the force 
level variances for the different initial force levels and different 
breakpoints for the case of symmetric parity. In Fig. 22 the variances 
for a battle with equal attrition rate coefficients but varying initial 
force levels show that H 4 is only true for the loser. This point is 
emphasized by Fig. 23, where with constant and equal initial force levels 
the attrition rate coefficients were varied. 
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Fifur* 20 - rORCE LEVEL VARIANCES OVER TIME 
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Vqrian 




JTi«ur« 21 - VARIANCES IN SYMMETRIC PARITY 
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Figure 22 - VARIANCES WITH DIFFERENT INITIAL FORCE LEVELS 
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Fi*ur« 23 - VARIANCES WITH DIFFERENT ATTRITION RATE 
COEFFICIENTS 
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V. COMPUTATIONAL ASPECTS 



A. USING THE NUMERICAL SOLUTION 

For most of the numerical work, the state probabilities have been 
obtained using the fourth-order RUNGE-KUTTA method. The accuracy of the 
results was increased by substituting the available analytical results 
for the Regions I, II and III. 

The disadvantage is that this solution method needs a lot of CPU-time. 
For the battle in Example 1, with the data given in Table 1, the calcula- 
tion of the state probabilities, expected force levels and variances led 
to a CPU-time on an IBM-360 computer of almost 90 minutes with a time 
step size of 0.05 minutes. That was the reason why, in the examples the 
attrition rate coefficients were increased by one magnitude, which brought 
the computer usage down to a CPU-time of around 12 minutes with the same 
time step size of 0.05 minutes. The sum of the state probabilities at 
every time step was used as a measure of accuracy. Surprisingly, its 
deviation from 1.0 was always in the fourth or fifth decimal, which proves 
the robustness of the RUNGE-KUTTA method. Also, the fact that for the 
Regions I, II and III analytical solutions were substituted at each time 
step did not change the final outcome considerably. It was found that 
without the analytical partial solutions the sum of probabilities tended 
to be slightly higher, changing the sum of probabilities from, for 
example, 0.99998 to 1.00003. 
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B. THE EQUAL ATTRITION RATE COEFFICIENT SOLUTION 

After the development of the equal attrition rate coefficient solution 
(EARCS) outlined in Ch. III.C, it was implemented on an IBM-360 comput- 
er. The CPU-time for the calculation of the state probabilities, expected 
force levels and variances for a given point in time always stayed below 
20 seconds, which emphasizes its computational advantage. 

Further research showed two shortcomings of the EARCS, which are 
easily overlooked. The coefficients C(m,n) vary over a wide range start- 
ing at 1.0 and depending on the initial force levels. 

For example for a battle with 



nig = 20 


m bp 


= 15 


O 

II 

o 

c 


n bp 


= 20 


f = 0.08 


t f 


= 1.60995 


C ( 20 ,40) = 1.0 


but 


C(16 ,21 ) = 4. 9297077276x1 0 33 



Also, the binomial coefficients in equations (3.23) and 3.24) vary over a 
wide range starting at 1.0 to 1352078 for the above example. This shows 
that the capacity the computer, in terms of the number of significant 
digits, is able to carry limits the implementation of this solution in 
the present form. 

Another reason why this solution is not the end of the numerical 
problems is the fact that in the given form, equations (3.23) and (3.24) 
require computation of an alternate sum consisting of terms whose factors 
are relatively large binomial coefficients and small numbers between zero 
and one. This produces truncation errors, which yield nonsensical results 
like negative variances and sums of probabilities greater than one. 
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VI. CONCLUSIONS 



Throughout this thesis, only LANCHESTER-type square-law attrition 
with fixed initial force levels and fixed attrition rate coefficients 
has been considered for a deterministic and a stochastic version. The 
stochastic version required much more computational effort. So, given 
the need for an analytical model as opposed to the use of simulation, 
there is not much to gain from the application of a stochastic model 
when the force levels are large and the forces are not near parity. In 
these cases the deterministic version essentially produces the same 
results, at least in qualitative terms. For smaller force levels or 
forces near parity, the stochastic version may be helpful to get more 
information about the dynamics of combat. 

Also, there exist one case where the deterministic version cannot 
differentiate between several different battles. That is the case of 
equal initial forces, equal attrition rate coefficients, but different 
breakpoint force levels. This also rectifies the further development of 
the equal-attrition-rate-coefficient-solution (EARCS). 

Using an analytical model, the discussed way of introducing random- 
ness into the model seems not to be very enlightening. Therefore, it is 
suggested that another way to include random effects should be explored. 
The author's opinion is that working with attrition rate coefficients 
which are random variables seems more promissing to gain insight into 
the dynamics of combat. Further down the line there should be some 
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consideration on the usage of combinations of the possibilities to 
include random effects as given in Chapter II. B. 
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16 
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APPENDIX A 



C THIS PROGRAM COMPUTES THE STATE PROBABILITIES USING THE FOURTH 
C ORDER RUNGE-KUTTA METHOD. PLOTS OF THE EXPECTED FORCE LEVELS AND 

C VARIANCES AS HELL AS 3-D-PLOTS FOR THE STATE PROBABILITIES ARE 

C OBTAINED. 

C 



DIMENSION PI (33, 33). P2 (33,33), T3D (20) 

DIMENSION XMT (500) , XNT (500) , OT (500) , VM (500) , VN (500) , 

1 VMB (551 ) , VNB (55 1 ) , VQ (551) . VARX(551), VARY (551), TIM (551) 
DIMENSION EST (313) , D (31 3) , F (2) , SI ZE (2) , KX (100) , KT (100) , HK (4 1 , 4 1 , 3) 
1 , PC (41 ,41) ,CST (41) ,CETERM (41) 



LOGICAL*! I DN (41,41) 

DIMENSION DETERM (313) 

REAL*8 TTL (12) /12*' '/ 

REAL K1 , K2, K3, K4 

CALL ERRSET (208. 600. - 1. 1 ) 

F (1) «0. 

F (2) *0. 

L I NES=»Q 
SIZE (1) *6. 

SIZE (2) -8. 

OFLAG • 0. 

NKXY* 1 00 

READ (5,21) MBP, MO, NBP, NO 

READ (5,22) AA. BB 

READ (5, 22) H, FT I M 

READ (5,21) N3D 

READ (5,22) EPSQ, EPSPTT 

T3D (1) - 9999. 

IF (N3D .LT. 1) GO TO 25 
READ (5,22) (T3D(I), I-1.N3D) 

25 CONTINUE 



C 

C N3D-NUMBER OF 3D-PLGTS, H=S1ZE OF TIME STEP 
C FTIME-FINALTIME, EPSQ, EPSPTT ARE ZERO LEVELS 

C T3D (I ) ■ FRACTION OF FT I ME WHEN TO PLOT 3D 

C 

DO 10 I-1.N3D 
T30 (I) *FT IM*T3D (I ) 

10 CONTINUE 

21 FORMAT (1615) 

22 FORMAT (8F10.5) 

WRITE (6,28) MBP. MO, NBP, NO 

28 FORMAT (// 5X, *M S 4 N S\ 5X. 416 /) 

WRITE (6,29) AA, BB 

29 FORMAT (/ 5X, 'A 4 B\ 9X. 2F10.3 /) 

WRITE (6,30) H, FT I M 

30 FORMAT (/ 5X, 'INITIAL H 4 FINAL TIME (« LOOPS * H) 5X, 2F 10.3/1 
WRITE (6,31) N3D 

31 FORMAT (// 5X, '* OF 3-D PLOTS'. 2X. 16 /) 

IF (N3D .LT. 1) GO TO 36 
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51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 
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66 

67 
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75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 



WRITE (6, 24) EPSQ, EPSPTT 

24 FORMAT!/ 5X. ’EPSQ* \FI0.5. 5X. ’EPSPTT- \ FlO. 5) 
WRITE (6, 32) (T30 ( I ) , 1 = 1. N3D) 

32 FORMAT!/ 5X. ’AT TIME’. 2X, 10F10.3 /) 

36 CONTINUE 
M01 = MO + 1 
N01 - NO ♦ 1 
NS - NO 1 - NBP 
MS « M01 - MBP 
MSL = MS + 1 
NSL - NS +1 
BMO » BB * MO 
ANO - AA « NO 
QANB = ANO / BB 
OBMA * BMO / AA 
8AT = - (BMO ♦ ANO) 

DO 47 Ml = 1 , MS 
VM (Ml) = MBP ♦Ml — 1 

47 CONTINUE 

DO 48 N1 = 1 , NS 
VN(N1) - NBP ♦ N1 - 1 

48 CONTINUE 
L - 0 

N3 » 1 

TPLOT - T3D (N3) 

TIME - 0. 

50 CONTINUE 

IF (TIME .GT. FT I M) GO TO 210 

52 CONTINUE 
TIME - H * L 
L » L + 1 

L8 - 0 
T1M(L) = TIME 
IF ( L . EQ. 2) GO TO 53 
IF (TIME .LT. TPLOT) GO TO 60 
C 

C SET NEXT PLOT TIME 
C 

N3 » N3 ♦ 1 
TPLOT = T30 (N3) 

IF (N3.GT.N3D) TPL0T=9999. 

53 CONTINUE 
C 

C L8 « L WHEN TIME TO PLOT 
C 

L8 - L 

E8T - EXP (BB * TIME) - 1. 

EAT = EXP (AA * TIME) - 1. 

AE - QANB * EBT 
BE » QBMA * EAT 
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1 22 

123 

124 

125 
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140 
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BAE - EXP (BAT * TIME) 

FJ - 1. 

JJ - 1 
C 

C NO VALUE OF J = ZERO AT MG, NO 
C 

FK - 1. 

KK - 1 

QFK-1. 

QF Je 1 . 

c 

C NO VALUE OF K = ZERO AT MO, NO 
C 

60 CONTINUE 

IF (L .EQ. 1) GO TO 500 
QQ=0. 

DO 400 N1 = 1 , NS 
N - N01 - N 1 
NL - NSL - N 1 
NL1 - NL ♦ 1 
ANJ - AA * N 
C 
C 

DO 300 Ml = 1 , MS 
M - M01 - Ml 
ML - MSL - Ml 



MLl 


a 


ML «■ 


1 








BMI 


at 


BB m 


M 








ABA 


a 


ANJ 


♦ BMI 








IF 


(M 


.EO. 


MO) 


GO 


TO 


502 


IF 


(M 


.EO. 


MBP) 


GO 


TO 


504 


IF 


(N 


.EO. 


NO) 


GO 


TO 


506 


IF 


(N 


.EO. 


NBP) 


GO 


TO 


507 



C 

C OEFAULTS TO ALL INSIDE POINTS 

C 

K 1 » AN J*P2 (ML 1 , NL) ♦ BMI *P2 (ML, NL 1 ) - ABA*P2 (ML . NL) 

PT - AN J*Q. 5* (PI (ML 1 , NL) +P2 (ML 1 , NL) ) ♦ BM1*0.5 * (PI (ML 
1 P2 (ML , NL 1 ) ) 

K2 * PT - ABA * (P2 (ML , NL) ♦ H*0.5*K1) 

K3 - PT - ABA * (P2 (ML , NL) ♦ H*0.5*K2) 

K4 - ANJ«P1CML1.NL) * BMI *P 1 (ML , NL 1 ) - ABA* CP2 (ML , NL ) ♦ 
Pi (ML , NL) * P2 1ML.NL) + CH/6.0) * (K1+2.0*K2 ♦ 2.0*K3 ♦ 
GO TO 200 

502 IF (N . EO. NO) GO TO 503 
IF (N .EO. NBP) GO TO 507 
C 

C M=MO, N^=NQ, NBP 

C 

K1 - BMI *P2 (ML , NL 1 ) - ABA*P2 (ML , NL) 



, NL 1 ) ♦ 



H*K3) 

K4) 
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151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 
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163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 



o o o 



PT - BMIxO.S x (PMML.NL1) + P2 (ML . NL 1 ) ) 

K2 - PT - ABA x (P2 (ML.NL) ♦ HkK1x0.5) 

K3 - PT - ABA x (P2 (ML , NL) + HxK2x0.5) 

KM =• BMI x Pl(ML.NLl) - ABA x (P2 (ML.NL) ♦ HxK3) 

P 1 (ML, NL) = P2 (ML , NL) + (H/6.) x (Kl+2.0xK2 + 2.0xK3 ♦ KM) 

IF (LB .EQ. L) GO TO 601 
GO TO 200 
C 

C M*MO, N=NO 

C 

503 CONTINUE 

PI (ML, NL) - EXPl-ABA x TIME) 

WRITE (6. J02) L. N. ML.NL. Pl(ML.NL) 

GO TO 200 

SOU IF (N .EQ. NBP) GO TO 505 
C 
C 

C M»MBP. N— • NBP 

c 

K1 - ANJ x P2 (ML 1 , NL) 

K2 ■ ANJ x 0.5 x (PI (ML 1 , NL) ♦ P2 (ML I . NL) I 
K3 - K2 

K4 - ANJ x PI (ML 1 , NL) 

PI (ML.NL) => P2 (ML.NL) ♦ (H/6.0) x CK 1 + 2.0xK2 ♦ 2.0xK3 ♦ KM) 
GO TO 200 
C 

C M=MBP, N=NBP 

c 

505 CONTINUE 
PI (ML.NL) » 0.0 
GO TO 200 

N-NO. M-=MBP, MO 

506 K1 - ANJxP2 (ML1.NL) - ABAxP2 (ML. NL) 

PT - ANJ x 0.5 x (Pi (ML 1 . NL) ♦ P2(ML1.NL)) 

K2 - PT - ABA x (P2 (ML.NL) ♦ HxO.SxKl) 

K3 - PT - ABA x (P2 (ML.NL) ♦ Hx0.5xK2) 

K« - ANJ x P 1 (ML 1 , NL) - ABA x (P2 (ML , NL) + HxK3) 

PI (ML.NL) = P2 (ML.NL) + (H/6.0) x CK 1 ♦ 2.0 xK 2 ♦ 2.0xK3 ♦ K4) 
IF (LB .EQ. L) GO TO 602 
GO TO 200 

C 

C N'NBP, 4 N=NBP. M=MO 

C 

507 CONTINUE 

HI - BMI x P21ML.NL1) 

K2 - BMI x 0.5 x (PI (ML , NL 1 ) + P2 (ML.NL 1 ) ) 

K3 = K2 

K1 - BMI x Pi (ML.NL1) 



76 



201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 
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219 
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222 

223 

224 
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226 
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228 

229 
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233 

234 

235 
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243 
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245 
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Pl(ML.NL) = P2 (ML * NL) + (H/6.0) * CK 1 ♦ 2.0*K2 ♦ 2.0*K3 ♦ K4) 
GO TO 200 
C 

C ALSO M=MO , N— =NBP, NO 

C 

601 CONTINUE 

IF MBE.LT. 1.) .AND. (KK.GT.9) ) GO TO 255 
C 

C COMPUTE PTT IF BE > 0 4 KK < 10 
C 

PTT-BAE 

00 250 1ND-1.KK 
XI NO* 1 NO 
PTT-PTTmBE/XIND 
250 CONTINUE 
GO TO 260 
255 PTT-O. 

260 CONTINUE 
KK * KK ♦ 1 

IF ( (A8S (PI (ML , NL) -PTT)) . LE . EPSPTT) GOTO 199 
GO TO 700 
C 

C M=MBP, N-NO 

C 

602 CONTINUE 

IF ( (AE.LT. 1.) .AND. (JJ.GT.9) ) GO TO 270 
C 

c 

C COMPUTE PTT IF AE > 0 4 JJ < 10 
C 

PTT-BAE 

DO 265 lND= 1 , J J 
X1ND-1ND 
PTT*PTT*AE/XIN0 
265 CONTINUE 
GO TO 275 
270 PTT«Q. 

275 CONTINUE 
JJ » JJ ♦ 1 

IF UABS (PI (ML.NL) -PTT) ). LE . EPSPTT) GO TO 199 

700 CONTINUE 
C 

C TIME REDUCED 1/2 
C 

L » L - 1 

WRITE (6,701) TIME, M. N. H 

701 FORMAT!// 5X, *H VALUE IS REDUCED BY HALF AT : TIME = \ F8.3, 

1 • M « \ 15, ’ N - \ 15, * FROM \ F8.3 //) 

H - 0,5 m H 

WRITE (6, 102) L.M.N.KK, PTT, P I (ML . NL) , QFK . QFJ.QQ.QFLAG 
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102 FORMRT (/ 2X, 41S, 7F12.5) 

QFLAG=0. 

GO TO 52 

199 PTT=P 1 (ML , NL ) 

200 CONTINUE 
QQ=QQ + Pl (ML , NL 1 

300 CONTINUE 

400 CONTINUE 
QQ1*ABS (QQ- 1 . ) 

IF (QQ1.LE.EPSQ) GO TO 401 
QFLAG* 1 . 

GO TO 700 

401 CONTINUE 
C 

C COMPUTE VRRX IT ) , VARY (T 1 , N3AR ( T I , MBAR (T) , Q (T) 
C 

SM8 = 0. 

SN8 = 0. 

SO * 0. 

SM2 * 0. 

SN2 * 0. 

DO 415 Ml = 1 * MS 

RM = VM (Ml 1 
RM2 * RM*RM 
DO 410 N 1 = 1,NS 
RN « VN (Nl) 

RN2 * RN*RN 
PT =* PI (Ml , Nl) 

P2 (Ml » Nl ) = PT 
SM8 = SM8 ♦ RM*PT 
SN8 a SN8 «* RNmPT 
SO = SO ♦ PT 
SM2 = SM2 ♦ RM2 * PT 

SN2 * SN2 «• RN2 m PT 

410 CONTINUE 
415 CONTINUE 

VM8 (L) = SMB 
VNB(L) * SNB 
VO CLI * SQ 

VRRX (L) * SM2 - SMBmSMB 
VART(L) = SN2 - SNB*SNB 
C 

C RETURN TO MAIN LOOP (50) IF NOT TIME TO PRINT 
C 

IF (N30 .LT. 1) GO TO 430 
IF (L8.NE.L) GO TO 430 
IF (L8.EQ.2) GO TO 430 
C 

C ADJUST X, Y, VECTORS FOR PLOT 
C 
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301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 
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319 
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321 
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330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 



00 15 K* 1 , MS 

15 EST (K) = K- 1 
00 16 K= 1 , NS 
OETERM (K) = K- 1 

16 CONTINUE 

00 115 1*1,41 
CST (I) -I 
CETERM (I ) * I 
00 115 J* 1,41 
PC(I,J)«O.G 

115 CONTINUE 
MP*4 1 

NP*4 1 
ALP* 15. 

BETA*30. 

PMAX*P1 (1,1) 

00 116 1*1, NS 

WRITE (6. 199) I , PI (1 , 1) .PI (1 , 1 ) 

199 FORMAT (’ \ 15, 10X.2F15.5) 

116 CONTINUE 

DO 18 1*1, MS 
DO 18 J* 1 , NS 

IF (PI (I , J) .GT.PMAX) PMAX=P 1 (I , J) 

18 CONTINUE 
C0NST=8. O/PMRX 
DO 19 1*1, MS 

N* I + MBP 
00 19 J* 1 , NS 
M= J+NBP 

PI (I, J) =CGNST*P1 (I, J) 

PC (N, M) *P 1 (I . J) 

19 CONTINUE 

CALL PLT301 (CST. MP. CETERM, NP. PC. RLP, BETR, F . TTL. SI ZE. WK. ION.KX.KT.N 
1KXT, LINES) 

C 

430 CONTINUE 
C 
C 

C PLOT HERE 

C RETURN TO MAIN LOOP 

C 

GO TO 50 

500 CONTINUE 

00 501 I * 1 , MS 
DO 501 J * 1 , NS 
PI (I, J) * 0.0 

501 CONTINUE 

PI (MS, NS) * 1.0 
GO TO 401 
C 
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C FINAL TIME REACHED 
C 

210 CONTINUE 

WRITE (6,211) TIME, FT I M, L 

211 FORMAT!// 5X, 'COMPUTED TIME ', F10.3. 5X, 'INPUT FINAL TIME ', 
1 F 10. 3, 5X, '• OF LOOPS TO REACH FINAL TIME ', 15 //) 

DO 215 I »1,L 

WRITE (6,213) VMB(I), VNB(l), VO ( I ) , VARX(I), VART(I), TIM(l) 

213 FORMAT (2X, 12F10.5) 

215 CONTINUE 

650 FORMAT ('1 ') 

WRITE (6,650) 

CALL PLOTP (TIM, VMB,L,0) 

WRITE (6,650) 

CALL PLOTP (TIM, VARX,L,0) 

WRITE (6,650) 

CALL PLOTP (TIM, VNB,L,0) 

WRITE (6,650) 

CALL PLOTP(TIM, VART,L,0) 

STOP 

DEBUG SUBCHK 
END 
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APPENDIX B 



C 

C THIS PROGRAM CALCULATES THE STATE PROBABILITIES FOR THE EQUAL 
C ATTRITION RATE COEFFICIENT SOLUTION (EARCS) . 

C A 3-D PLOT IS PRODUCED USING THE VERSATEC PLOTTER. 

C 

IMPLICIT REAL*8 (A-H.O-Z) 

CALL ERRSET (208.256. 10. I) 

DIMENSION C (50.501 .M(50) .N (50) , F (50) , PTMN (50.50) 

REAL*4 SIZE (2) . FL (2) .WK (41.41. 3) . X (41 ) . T (41) .P (41.41) 
DIMENSION KX (100) .KT (100) 

LOGICALhI 1DN (41 .41) 

REAL*8 TTL (12) /12«* V 

READ (5. 100) MO.MBP.NO.NBP 
READ (5. 101) A 
READ (5. 102) TIME 
READ (5.103) (F (I) . 1 = 1,41) 

C 

C F (I) * 1 /I-FACTQR1 AL , DONE TO SPEED UP THE PROGRAM 
C 

WRITE (6,802) MO. MBP. NO . NBP, A , T I ME 
802 FORMAT ( f * , 4 1 5. 2F 1 0. 5) 

MD=M0-MBP 
ND=NO-NBP 
M01=M0*1 
N01 = N0M 
MBP1=MBP*1 
NBP1=NBPM 
MBP2=MBP ,f 2 
NBP2=NBP*2 
RMO=M0 
RN0=NQ 
DO 10 1=1,41 
DO 10 J= 1 , 4 1 
C Cl, Ji =0.0 
PTMN (I, J) =0.0 
PCI, J) =0.0 
10 CONTINUE 

C CMOl.NOl) =1.0 
DO 20 1=2, MD 
J=MBP+I 

C C J, NO 1 ) = RNG** (MD- 1 ♦ 1 ) 

20 CONTINUE 

DO 21 1=2, ND 
J=NBP*1 

C CM01, J) * RMO** (ND- I ♦ 1 ) 

21 CONTINUE 

DO 22 1=2, MD 
MC=M01 - I 
MM=M01 - I ♦ 1 
MPLUS=MM^1 



81 



51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

66 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 



DO 22 J=2, NO 
NC=N01 - J 
NN*N01-JM 
NPlUS=NN+ 1 

C (MM, NN) *NC* C (MPLUS , NN) *MC* C (MM.NPLUS) 
22 CONTINUE 
C 

C END OF COEFFICIENT CALCULATION 
C 

F1=-A*TIME 
DO 25 1 = MBP2, MO 1 
DO 25 J=NBP2, NO 1 
1 J=I ♦ J 

K=M0+NO- I J *2 

F2* (1 .O-DEXP (FI)) **K 

KK=K+1 

I J=1 J-2 

F 3=DEXP (F1*I J) 

PTMN (I , J) =F2*F3*C (I , J) *F (KK) 

25 CONTINUE 

100 FORMAT (415) 

101 FORMAT (FlO. 5) 

102 FORMAT CF 1 0 . 5) 

103 FORMAT CD 1 7 . 1 I) 

C 

C FOR NBP<N<NO 
C 

DO 30 11=1, ND 

NF 0RCE=NBP+ 1 1 

J*M0+N0-MBP-1-NF0RCE 

FACT=MBP+1. O+NFORCE 

SUM1 = (1 . O-DEXP (F 1 *F ACT) ) / ( A*FACT) 

SUMI NT=$UM1 
DO 301 K* 1 , J 
F ACT =FACT* 1 . 0 
ADDFAC* (-1.0) k*K 

FRACTN* (1. O-DEXP (F 1 *FACT) ) / ( A ♦'FACT) 

COMBT = 1.0 
DO 302 KJ= 1 , K 
RKJ*KJ 

COMBI = (J-RKJ+1 .0) /RKJ 
COMBT* COMBT x COMB I 
302 CONTINUE 

SUMINT=SUMINT*ADDFAC*FRACTN*COMBT 
301 CONTINUE 
NMBP= J + 1 
NB0UND=NFORCE*1 

PFAC=A*NF0RCE*F (NMBP) *C (MBP2, NB0UND) 

PTMN (MBP 1 , NBOUND) =PFAC*SUMINT 
30 CONTINUE 
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101 

102 

103 

104 

105 

106 

10 ? 

108 

109 

110 

111 

112 

113 

114 

115 

116 
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123 

12*1 
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134 

135 

136 

13? 

138 

139 

mo 

141 

142 

143 

144 

145 

146 

147 

146 

149 

150 



CJ (J 



c 

C FOR MBP<M<MO 
C 

DO 40 11=1. MO 

MFQRCE-MBP+Il 

J»M0*N0-NBP-1 -MFORCE 

FACT-NBP+l .O+MFORCE 

SUM1- (1 .0-OEXP (F 1 wFflCT) i / ( A-FACT) 

SUMlNT-SUMl 

00 401 K= 1 , J 

FACT-FACT+l .0 

ADDFAC- (-1.0) *«K 

FRRCTN= (1.0-DEXP (F 1 *FACT) ) / ( A*FACT) 
C0MB7=1 .0 
DO 402 K J- 1 , K 
RKJ=K J 

C0MBI= < J-RKJ+ 1 . 0) /RKJ 
COMBT-COMBT *C0MB1 
402 CONTINUE 

SUH1NT=SUM1NT + ADDFAC * FR ACT NwCQMBT 
401 CONTINUE 
MNBP- J + 1 
MBOUND-MFORCE+1 

PFAC*A*MFORCE*F (MNBP) *C (MBOUND, NBP2) 

PTMN (MBOUND. NBP1) =PFAC«SUHINT 
40 CONTINUE 

DO 45 1 =MBP1 . M01 
DO 45 J=NBP1 , N01 
K=I-1 
L-J-l 

WRITE (6,801) K.L, PTMN (I. J) 

801 FORMAT ( * PIT.*. 15, \ M5. ’) = *. D17.ll) 
45 CONTINUE 

&ATA ADJUSTMENT FOR PLOT 

NR0W=41 
NC0L=4 1 
NKXT=100 
L1NES=0 
ALPHA-15. 

BETA-30. 

FL (1) -0.0 
FL (2) -0.0 
SIZE (1) =6.0 
SIZE (2) =8.0 
C 

C SCALING 
C 

PMAX-PTMN (1 . 1) 
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188 
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197 
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199 

200 



50 



DO 50 1*1,41 
X (I) *1 

DO 50 J= 1 , 4 1 
T (J) *J 

I F (PTMN ( I , JJ .GT.PMAX) PMAX=PT MN ( I , J) 

CONTINUE 
C0NST*8 . Q/PMAX 
DO 51 1*1,41 
DO 51 J* 1 , 4 1 

P(I,J) * SNGL (CGNST*PTMN (I, J) ) 

51 CONTINUE 
C 
C 

C IF OTHER FORCE LEVELS CHANGE PLOT ARGUMENTS 
C 

CALL PLT3D1 (X, NRGW, T, NCOL , P, ALPHA, BETA, FL, TTL, SI ZE, WK, 

* ION, KX.KT,NKXT, LINES) 

C 

C EXPECTED VALUES AND VARIANCES 
C 

EM*0. 0 
EMM=Q. 0 
VARM*0. 0 
SPR0B*0. 0 
DO 60 1 I *MBPl , MOl 
1*11-1 
SPT*0. 0 

DO 601 J J=NBP1 , NO 1 
SPT*SPT ♦ PTMN (I I , JJ) 

601 CONTINUE 
EM=EM+!*SPT 
EMM*EMM+ I*I*SPT 
SPRGB=$PRGB+SPT 
60 CONTINUE 

VARM*EMM-EM*EM 
EN*0.0 
ENN*0. 0 
VARN=0. 0 

DO 70 1 1 *NBP1 , N01 
1 * 11-1 
SPT*0. 0 

DO 701 J J = MBPl , MO 1 
SPT=SPT+ PTMN (JJ, II) 

701 CONTINUE 
EN=EN+I*SPT 
ENN=ENN+ I*I*SPT 
70 CONTINUE 

VARN*ENN-EN«EN 

WRITE (6, B05) EM, EN , VARM. VARN 

805 FORMAT ( * \ ’EXPECTED VALUES M , N\ 50X, * VARIANCE M.N ’,/ 
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201 

202 

203 

20M 

205 



h* \DI7.II.10X. D17.11.T51.D17. 11. 10X. Dl7.ll) 
HR1TE (6. 1100) SPROB 

1100 FORMAT ( ' SUM OF PROBABILITIES \D17.I1) 

STOP 

END 
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